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SOME  BUBBLE  AND  CONTACT  PROBLEMS 

JOSEPH  B.  KENT  R 

Abstract.  A  number  of  problems  invoking  bubble*,  in  a  fluid.  or  contact  of  surfaces.  or  both  arc 
considered  In  each  case  the  size  and  or  the  shape  ot  the  bubble.  or  the  location  ot  the  points  of  contact,  are 
unknown  in  advance  and  must  be  found.  When  the  problems  ate  tot u>ulaud  mathcnuticalh.  bounda;\ 
conditions  must  be  satisfied  at  the  bubble  surface  and  at  the  contact  points,  which  are  thus  free "  boundaries 
The  problems  are  post  buck  ling  behavior  of  an  elastic  tube  contact  problems  involving  a  buckled  elastka, 
steep  capillary  waves  with  trapped  bubbles,  deformation  of  a  bubble  in  a  uniform  flow,  distortion  of  a  bubble 
in  a  straining  flow  ,  free  oscillation  of  an  underwater  explosion  bubble,  and  forced  oscillation  of  a  bubble  in  a 
sound  field. 


1.  Introduction.  Physical  problems  involving  bubbles  or  drops  in  a  fluid  generally 
require  the  determination  of  the  size  and  or  the  shape  of  the  bubble  or  drop  surface. 
Problems  involving  the  contact  of  elastic  or  fluid  surfaces  usually  require  the  location  of 
the  points  of  contact.  In  both  cases,  boundary  conditions  must  be  satisfied  at  the 
unknown  bubble  or  drop  surface  or  contact  point.  Thus  these  surfaces  or  points  are 
“free"  boundaries,  which  must  be  found  as  part  of  the  solution  of  the  problem.  In  the 
case  of  bubbles  in  contact  with  themselves  or  with  other  surfaces,  there  are  contact 
points  on  the  bubble  surface.  These  points  are  thus  free  boundaries  of  free  boundaries, 
which  might  be  called  "free  free  boundaries",' 

W  ith  the  help  of  a  number  of  coworkers,  I  have  analyzed  various  free  boundary 
problems  involving  bubbles  and 'or  contact  of  surfaces.  1  shall  describe  some  of  them  in 
this  report,  emphasizing  the  reasons  for  considering  them,  the  principles  involved  in 
their  formulation,  and  the  nature  of  the  results.  In  all  cases,  some  analytical  results  and 
some  numerical  results  were  obtained.  I  shall  present  some  of  the  numerical  results  and 
some  figures  based  on  them,  indicating  what  the  numerical  methods  were  when  that 
seems  appropriate.  I  shall  just  mention  the  analytical  results,  referring  for  their  detailed 
presentation  to  the  relevant  published  articles. 

The  goal  of  this  report  is  to  show  that  certain  problems  involv  ing  contact  ot  surfaces 
can  be  formulated  and  analyzed  in  a  relatively  simple  way  In  some  cases,  the 
maintenance  of  contact  is  required  to  prevent  the  unphysical  phenomenon  of  a  surface 
crossing  over  itself.  When  the  surface  is  a  fluid-fluid  interface,  the  resulting  contact  mav 
lead  to  the  occurrence  of  bubbles  or  drops  When  the  bubbles  or  drops  are  deformed, 
parts  of  their  surfaces  may  again  come  into  contact  with  one  another  and  form  new 
bubbles  or  drops.  Since  phenomena  involving  the  contact  of  surfaces  are  of  widespread 
occurrence,  it  is  to  be  hoped  that  the  results  reported  here  will  stimulate  further 
investigations  of  them. 

2.  Postbuckling  behavior  of  an  elastic  lube.  While  investigating  blood  flow  m 
arteries  and  veins,  mv  colleague  S.  I.  Rubinow  and  I  realized  that  the  pressure  of  the 
blood  in  a  vein  could  fall  below  the  ambient  pressure  outside  the  vein.  As  a 
consequence,  the  wall  of  the  vein  could  buckle,  thus  changing  the  vein  cross  section 
from  circular  to  some  other  shape,  and  thereby  increasing  the  vein  s  resistance  to  blood 
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flow.  Therefore,  with  J.  E.  Flaherty  [2],  we  studied  the  buckling  and  postbuckling 
behavior  of  an  elastic  tube. 

We  assumed  that,  for  any  given  pressure  difference  p,  the  equilibrium  configura¬ 
tion  of  the  tube  would  be  a  cylinder,  and  we  denoted  by  C'  its  cross  sectional  curve.  Then 
the  equations  of  elastostatics  governing  the  tube  are  the  same  as  those  of  an  elastic  ring 
C  compressed  by  a  force  p  per  unii  length.  This  problem  had  been  studied  by  I. 
Tadjbakhsh  and  F.  Odeh  [3].  They  assumed  that  the  ring  is  inextensible  and  is  governed 
by  the  Euler-Bernoulli  beam  theory,  with  the  bending  moment  proportional  to 
Jfcl.v  i  -  1.  Here  kis  <  is  the  curvature  of  C  at  arclength  v  along  it.  and  1  is  the  curvature 
when  C  is  a  circle.  In  dimensionless  variables  these  assumptions  lead  to  the  equation 

(2. Ip  k"  +  \k  *  -  ck  -p  =  II,  0gs227r. 

Here  c  is  a  constant  to  be  determined.  In  addition,  k  must  be  periodic  with  period  2n. 
and  the  integral  of  k  over  a  period  must  be  2 77.  This  latter  condition  follows  from  the 
fact  that  kis)  =  rt  i.si,  where  fh.vi  is  the  angle  between  the  tangent  to  C  and  some  fixed 
direction. 

For  every  value  of  p  *  0.  this  problem  has  a  solution  k  =  1  and  c  =  \-p. 
corresponding  to  C  being  a  circle.  As  Tadjbakhsh  and  Odeh  [3]  showed,  this  unbuckled 
solution  is  the  only  solution  for  p  •£  3,  while  for  p  3  there  are  other  solutions.  Each 
other  solution  is  periodic  with  period  2 7 r/n  for  some  integer  n  2  2  and,  for  a  suitable 
origin  of  v.  it  is  odd  in  v.  Thus  each  such  solution  satisfies  the  conditions 

1 2.2 1  k  1(11=0,  k  1 77  >1 1  -■ ■  0.  k  1  v  1  ds  -  it ■  n. 

Furthermore,  it  exists  onk  for  p  ■  p  „.  where  p ...  is  called  the  mh  buckling  load  They 
computed  the  solutions  numertcallv  for  the  first  few  values  of  n  for  .1  range  of  p  *p„„. 
For  the  nth  solution,  when  p  reached  a  second  critical  value  p  ...  »t  2.  n  pairs  of  points 
came  into  contact  with  one  another,  while  for  n  =  2.  one  pair  of  opposite  points  came 
into  contact  For  p  ■  p  „.  the  ring  crossed  over  itself. 

Crossing  over  may  be  possible  for  opposite  parts  of  a  ring,  if  one  part  comes  out  of 
the  plane  containing  the  rest  of  the  ring  However  it  is  not  possible  for  a  tube  Therefore 
we  reexamined  the  problem  for  the  nth  solution  with  p  •  p  „.  to  find  a  solution  which 
does  not  cross  over  itself  W  e  required  that  for  p  ■  p  one  pair  of  points  in  each  period 
be  in  contact  with  one  another  For  n  2.  these  points  are  the  endpoints  of  the  period, 
so  there  is  onlv  one  pair  in  contact,  and  thev  are  at  >  in  2  For  n  2  the  locations  ±x; 
of  the  contact  points  must  be  determined 

In  the  case  n  2  we  found  that  the  second  condition  in  >2  2 1  had  to  be  changed  so 
that  2  2  '  became 


2  1  i  ill  o  i  -  2 


( 


k  •  \  *  il\  •  7t  2  for  p  •?  p  ■. 


We  solved  2  I  and  2  2  numen.  allv  to  obtain  k  >  v  >  and  <  for  v arious  v alues  of  p  from 
p  p.  *  up  to  ;  p  '  24 '  I  hen  we  solved  2  1  1  and  <  2  * 1  for  p  from  p  ~  p  7  = 
s  24’uptop  p  I"  <4  when  the  curvature  k  ■  n  2 1  at  the  contact  point  became  zero. 
F01  p  p  *c  assumed  'hat  the  solution  has  an  interval  of  contact  extending  from 
wmoilui  •  "  up  to  ~  2  I  hen  we  showed  that  v  and  the  shape  of  the  remaining  part 

of  the  v-luti  m  .ouid  be  obtained  bv  .1  similarity  transformation  from  the  result  for 
p  p  f  he  results  are  shown  in  Fig  1 .  and  corresponding  results  for  n  =  4  are  shown  in 
F'g  2 
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Once  we  had  the  shape  of  the  tube's  cross  section  for  n  =  2  and  each  value  of  p.  we 
solved  the  equation  An  =  -1  within  the  cross  section  with  u  =  0  on  the  boundary  for 
each  p.  This  yielded  the  velocity  of  the  flow  of  a  viscous  fluid  along  the  tube.  The  results 
showed  that  the  conductivity  of  the  tube,  i.e.,  the  flux  of  fluid  for  unit  pressure  gradient 
along  the  tube,  decreased  drastically  when  the  tube  began  to  buckle.  These  results  have 
been  confirmed  quantitatively  by  the  experiments  of  H.  Zeller  and  R.  Wirtz  [4],  They 
measured  the  flow  through  a  rubber  tube  buckled  by  a  pressure  difference  p,  and  their 
results  agreed  fairly  well  with  our  calculations,  as  is  shown  in  Fig.  3.  Their  correspond¬ 
ing  results  for  a  tube  of  polyvinylchloride  (PVC),  also  shown  in  Fig.  3,  do  not  agree  so 
well,  although  they  would  do  so  if  the  Young's  modulus  E  used  in  scaling  the  horizontal 
axis  were  replaced  by  E/2.  The  cross  sections  of  buckled  Penrose  tubing  in  various 
stages  of  collapse  during  flow  look  like  the  cross  sections  in  Fig.  2  (see  [5],  [6]),  as  do  the 
cross  sections  of  the  vena  cava  in  normal  breathing  [7], 


Flo.  3.  The  normalized  flow  through  a  tube  as  a  function  of  the  pressure  difference  between  the  outside  and  the 
inside.  The  open  circles  and  the  dashed  curve  are  based  on  the  theory  described  in  §  2.  The  black  circles  are 
experimental  values  for  a  rubber  tube  measured  by  H.  Zeller  and  R.  Wirtz  [4].  The  black  triangles  are  their 
measured  values  fora  tube  made  of  polyvinylchloride  {PVC).  They  would  have  agreed  better  with  the  theory if  a 
smaller  value  had  been  used  for  Young's  modulus  E  of  PVC. 


3.  Contact  problems  (or  a  buckled  elastica.  In  1744,  Euler  formulated  and  solved 
the  problem  of  the  buckling  of  a  thin  elastic  rod  or  column,  which  he  called  an 
"elastica",  subject  to  a  compressive  load  P  at  its  ends.  He  found  that  a  rod  clamped  at  its 
ends  remained  straight  or  unbuckled  for  OS  PsPbi  =  4  w2,  where  Phi  is  the  dimension¬ 
less  first  buckling  load.  For  P  >  Pb\  the  rod  could  buckle,  and  as  P  increased  it  deformed 
more  and  more  until  a  pair  of  points  came  into  contact  at  some  value  P.  i .  For  P  >  P.  i  the 
elastica  crossed  over  itself.  If  we  consider  the  elastica  to  lie  in  a  plane,  or  to  be  a  cross 
section  of  a  cylindrical  elastic  sheet,  it  cannot  cross  itself.  Therefore  J.  E.  Flaherty  and  I 
sought  a  solution  for  P>  Pt  ]  which  did  not  do  so  [8]. 

In  dimensionless  variables,  we  let  0(s)  be  the  angle  between  the  x-axis  and  the 
tangent  to  the  elastica  at  arclength  position  s,  with  s  ranging  from  0  to  1 .  Then  the 
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equation  of  equilibrium,  according  to  the  Euler-Bernoulli  beam  theory,  is 

(3.1)  0„+Psin0  =  O,  Ogygl. 

The  condition  that  the  ends  are  clamped  along  the  .t-axis  yields 

(3.2)  0(0)  =  0,  0(  1 )  =  0,  [  sin  0(s)  rfs  =  0. 

Jo 

We  solved  (3.1)  and  (3.2)  numerically  for  various  values  of  P  from  Phl=4n 2  to 
P.  i  =72.187,  at  which  contact  first  occurs.  (See  Fig.  4.i 


0  00  .25  50  75  1.00 


P*  40  000 


P*  200  000  P*400  000  P=600  000 

*,=0.098  *,=0.076  *,  =  0.065 


Fig  4.  The  shape  of  a  clamped  elastica  in  the  mode  n  =  1  for  six  values  of  the  end  load  P.  Contact  occurs  at 
P,  =  72  187  For  P  >  P£  the  similarity  solution  was  used  to  compute  the  shape  of  the  loop  between  the  points  of 
contact. 
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For  P>  P.  i,  we  require  contact  to  occur  between  the  two  points  s  i  and  1  -  st,  with 
s,  to  be  found  in  the  range  0  <  Si  <  i.  Then  (3.1)  holds  only  in  the  intervals  0  <  s  <  t,  and 
1  -  v,  <  s  <  1 ,  while  between  these  intervals  it  is  replaced  by 

(3.3)  +  (P- R)  sin  ft  ~  0,  s,<s<l-si. 

Here  R  is  the  contact  force,  which  is  to  be  found.  We  seek  a  solution  in  which  0(s)  is 
antisymmetric  about  s  =  ],  so  it  suffices  to  solve  (3.1)  for  0<s<s,  and  (3.3)  for 
Si<S<\.  In  addition,  we  have  the  six  conditions 

0(O)  =  O,  0(Si)  =  tt/2,  0(\)-O,  I  cos0(s)ds=O,  6  and  0,  continuous  at  Si . 

(3.4) 

These  six  conditions  suffice  to  determine  S\,R  and  the  two  integration  constants  in  each 
interval. 

We  found  that  the  solution  of  (3.3),  which  describes  the  closed  loop,  could  be 
determined  for  any  P>PC\  from  the  solution  at  Pc  \  by  a  similarity  transformation. 
Therefore  for  P>P,  ,  we  only  had  to  solve  (3.1)  numerically,  together  with  the 
appropriate  conditions  from  (3.4). 

The  results  for  three  such  values  are  shown  in  Fig.  4.  We  also  found  that,  for  P 
large,  5,  is  given  by 

(3.5)  5,-  1.854/’  '  --6.055P  P»  1. 

Thus  the  contact  points  approach  the  endpoints  as  P  tends  to  infinity. 

These  results  can  be  confirmed  qualitatively  with  a  strip  of  flexible  plastic,  by 
holding  one  end  in  each  hand  and  bringing  the  ends  toward  one  another. 

4.  Capillary  waves  with  trapped  bubbles.  Sufficiently  short  waves  on  the  surface  of 
a  liquid  are  governed  mainly  by  surface  tension.  For  many  years  such  capillary  waves 
were  studied  on  the  assumption  that  the  amplitude  a  was  much  less  than  the  wavelength 
A,  so  that  the  steepness  s=a/\  was  small.  Then  it  was  possible  to  linearize  the 
governing  equations.  In  1957,  D.  G.  Crapper  [9]  sought  to  express  the  wave  motion  as  a 
power  series  in  s.  From  the  regularity  of  the  first  few  coefficients  he  guessed  the  general 
term,  verified  it,  and  then  was  able  to  sum  the  series  in  closed  form.  His  result  yielded 
waves  ranging  from  sinusoidal  waves  for  5  small  up  to  a  steepest  one  at  s  =  .73,  when 
adjacent  waves  touched  one  another  at  one  point.  For  5  >  .73,  his  result  was  unphysical 
since  adjacent  waves  overlapped  one  another. 

In  view  of  the  studies  reported  in  the  two  preceding  sections,  I  believed  that  it 
should  be  possible  to  obtain  steeper  waves  by  enforcing  contact  of  adjacent  waves  and 
preventing  overlap.  Therefore  J.-M.  Vanden-Broeck  and  I  [10]  did  just  that.  Following 
Crapper,  we  considered  steady  two-dimensional  periodic  capillary  waves  on  the  surface 
of  an  inviscid  incompressible  fluid  of  infinite  depth.  We  considered  the  complex  position 
x  +  iy  as  a  function  of  the  complex  potential  ip+itb  in  the  halfplane  i^SO.  The 
streamline  ib  =  0  corresponded  to  the  free  surface.  At  it,  the  pressure  in  the  fluid  had  to 
exceed  the  pressure  in  the  air  above  the  surface  by  the  surface  tension  times  the 
curvature  of  the  surface.  In  dimensionless  variables  this  led  to  the  equation 

(4.1)  (xl  +  >•*)”'-  1  -  T(xvy„  +  y*rVJ  =  0  on</f  =  0. 


Here  T  is  a  dimensionless  constant  proportional  to  the  surface  tension.  The 
specification  of  the  fluid  velocity  at  y  =  -oo  is  that  (xtf,  y„)  -*  (1, 0)  as  ib  -*  -oo. 


To  complete  the  formulation  of  the  problem,  we  employed  the  assumed  periodicity 
in  v  of  the  solution  and  the  Hilbert  transform  to  obtain 


v  i c. lb 


f  [a.  i  v'  .  *b  1  Ifcot  ~i  c 


cos  ->  e  -  v- 1 ;  Jc 


In  terms  of  v,,  the  steepness  v  is  given  by 


v.  i  v  •  H'de 


Equations  t4.lt-t4.3i  constitute  a  problem  for  the  determination  of  the  periodic 
functions  vv.(y,  0).  v.-ty,  Ot  and  T  for  a  given  value  of  v. 

When  j  is  large  enough,  we  modify  the  preceding  formulation  by  requiring  contact 
to  occur  at  the  point  y  =  a.  i h  =  0  on  the  free  surface,  with  a  to  be  found  in  the  interval 
0 <  a  <  !.  Then  at  the  contact  point  we  require 


.v -to.  0)  =  0. 


Aim,  til 


Since  contact  of  adjacent  waves  results  in  a  trapped  bubble  of  air.  we  permit  the 
pressure  in  the  bubble  to  exceed  the  pressure  in  the  air  above  the  surface  by  the 
..mensionless  amount  P.  Then  we  must  replace  the  right  side  of  \4.1  >  bv  P  in  the  ranee 
a<*s\. 

We  solved  the  problem  without  contact  by  introducing  V  -  1  equally  spaced  mesh 
points  y.  in  the  interval  0  S  y  g  1 ,  the  v  alues  a  ,’  =  a.  ly ..  lb.  and  the  values  v  of  v.  '  y.  <b 
midway  between  adjacent  mesh  points.  By  using  the  trapezoidal  rule  in  4.2>  we 
expressed  the  y’,  in  terms  of  the  v  '.  Then  we  used  centered  four-point  interpolation  and 
difference  formulas  to  evaluate  the  derivatives  in  4.1  .  In  this  wav  we  converted  4  1 
into  a  set  of  X  -  1  nonlinear  algebraic  equations  for  the  \ and  i4.3unto  an  V  ~2md 
equation.  We  solved  these  equations  by  Newton  s  method  for  T  and  the  V  1 
quantities  v  using  ,V  =  40  and  .V  bo.  The  two  sets  of  results  agreed  to  four  decimals 
with  each  other  and  with  Crupper  s  solution. 

To  find  solutions  with  contact,  we  introduced  the  two  extra  unknowns  a  and  P.  and 
the  two  extra  equations  (4.4 1.  We  also  forced  n  to  be  a  mesh  point  by  separately 
subdividing  the  intervals  0  y:  y  %  «  and  4SyS;,  Then  we  solved  the  resulting  equa¬ 
tions  for  <>.  P.  T  and  the  r!  by  New  ton's  method  using  60  mesh  points.  For  each  v  alue  of 
v  we  used,  as  the  initial  approximation  to  the  solution,  the  solution  for  a  slightly  smaller 
value  of  s,  starting  with  Crapper's  solution  at  r  =  0.730. 

We  were  able  to  find  solutions  with  v  •  0.730.  and  we  computed  them  up  to  v  =  1 .5. 
It  seems  clear  that  we  could  compute  them  for  any  larger  value  of  s  by  increasing  the 
number  of  mesh  points.  Typical  solutions  are  shown  in  Fig.  5.  We  also  found  solutions 
with  s<  0.730,  in  the  interval  0.663  <  s  <  0.730.  so  that  two  solutions  exist  in  this 
range — Crapper's  and  this  new  one.  As  v  decreases  to  0.663.  the  curvature  of  the 
surface  at  the  contact  point  decreases  to  zero,  which  is  why  this  family  ends  there.  On 
the  other  hand,  as  s  -*  x,  each  bubble  becomes  infinitely  long  and  narrow,  with  its  w  idth 
proportional  to  e  Simultaneously  the  surface  above  them  tends  to  a  scries  of 
semicircles  tangent  to  one  another. 

5.  Deformation  of  a  bubble  in  a  uniform  Row.  Surface  lension  forces  a  bubble  in  a 
liquid  at  rest  to  become  spherical  in  three  dimensions  and  circular  in  two  dimensions. 
But  pressure  forces  due  to  a  fluid  flowing  around  a  bubble  can  deform  it  to  another 
shape.  We  have  seen  an  example  of  this  in  the  previous  section,  where  the  bubbles  were 
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long  and  tapered  at  their  upper  ends.  Hut  in  that  ease  each  bubble  was  affected  bv  the 
presence  of  the  upper  surface  and  the  other  bubbles.  It  is  more  important  to  know 
how  an  isolated  bubble  is  deformed,  since  bubbles  often  occur  in  isolation.  Therefore 
J.-M.  Vanden-Broeek  and  I  [  1 1  ]  investigated  this  deformation. 

We  considered  a  two-dimensional  bubble  in  a  steady  potential  flow  of  an  in  viscid 
incompressible  fluid,  ihe  flow  was  assumed  to  be  umtoiin  at  infinity  and  the  pressure 
within  the  bubble  was  assumed  to  be  constant.  On  the  unknow  n  bubble  boundary,  the 
fluid  velocity  was  required  to  be  dnecied  tangentially  and  the  bubble  pressure  was 
required  to  exceed  the  fluid  pressure  by  the  surtaee  tension  multiplied  by  the  curvature 
of  the  boundary,  upon  introducing  dimensionless  variables,  we  found  that  the 
coefficient  of  surface  tension  o  could  be  eliminated  by  using  hr/pU'  as  the  unit  ol 
length,  where  p  is  the  density  ol  the  fluid  and  U  is  the  velocity  at  infinity.  Then  the 
problem  contained  a  single  dimensionless  parameter  y.  defined  by 

(5.1)  y  =  2l  pi, -p,  \/pU2. 


Here  ph  is  the  pressure  in  the  bubble  and  p,  is  the  stagnation  pressure,  which  equals  the 
pressure  at  infinity  minus  pU'i 2. 

Before  solving  the  problem,  we  required  that  the  flow  and  the  bubble  be  symmetric 
about  the  r-axis.  which  we  took  parallel  to  the  velocity  at  infinity.  Then  it  sufficed  to 
consider  the  flow  outside  the  bubble  in  the  halfplane  >•«£().  This  flow  region  cor¬ 
responded  to  the  halfplane  i b  £  0  of  the  complex  potential  plane,  and  the  bubble  surface 
corresponded  to  the  segment  --  b  <  y  <  b  of  the  axis  y  =  0,  with  b  unknown.  It  was 
convenient  to  set  y  -  by  so  that  the  bubble  surface  corresponded  to  the  segment 
-l<y<l,  ib-  0.  Then  we  sought  i  -  iv  as  an  analytic  function  of  hy  +  f<A  in  this 
halfplane.  The  velocity  at  infinity  was  U„.  y, )  =  in.  Ol. 

By  applying  Cauchy's  theorem  to  v  -  iv  -  b  and.  taking  the  real  part,  we  obtained 
the  Hilbert  transform  relation 

1  ('  i  (y,0l 

(5.2)  .v,.ty,  Ol  =  />+-  I  — : - <ly  ,  -  1  <  y  <  I. 

n  J  ]  *£ 

We  also  wrote  the  pressure  jump  condition  on  the  bubble  surface  in  the  form 

(5.3)  ‘  =(>vvw  -A^y^X.r;.  '-y«l  '  ‘  -  y.  -l<yvl.  t//  =  <>. 

In  addition,  the  symmetry  about  the  r-axis  yielded 


(5.4) 


lim 


>T(y,  0) 


•‘i  x^ty,  0) 


=  oo 


These  equations  (5.2)— (5.4)  constituted  the  problem  for  the  determination  of 
xyty,  0),  \\  (ip,  0)  and  b  for  a  given  value  of  the  parameter  y. 

To  solve  these  equations,  we  required  the  bubble  to  be  symmetric  about  the  v-axis 
and  therefore  about  the  y-axis.  Then  we  could  restrict  our  considerations  to  the  interval 
OSy  <1  with  >y(0,  0)  =  0.  Within  this  interval  we  introduced  the  new  variable  a  in 
place  of  y  by  the  equation  y  =  1  -  a " .  This  eliminated  from  the  solution  the  singularities 
at  the  stagnation  point  y  =  1,  <b  =  0.  Once  this  new  variable  was  rC.med,  we  introduced 
mesh  points  and  values  of  the  unknown  functions  at  the  mrsh  points  just  as  in  the 
previous  section.  By  proceeding  in  the  manner  described  there,  we  converted  the 
equations  above  into  a  set  of  nonlinear  algebraic  equations  for  b  and  the  values  of  vy  at 
the  mesh  points. 
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To  sohe  these  equations  hv  New  ton  'method  for  a  ;n en  value  of  v.  we  needed  an 
initial  approximation  to  the  solution.  Therefore  we  determined  the  asymptotic  form  of 
the  solution  for  y  large,  when  the  bubble  tends  to  a  circle  of  radius  y  We  used  this  as 
the  initial  approximation  foi  a  Urge  value  of  s  and  iterated  until  we  obtained  a  solution. 
Then  ve  used  that  solution  as  the  initial  approximation  for  a  slightly  smaller  v  alue  of  y. 
and  so  on. 

Three  of  the  resulting  bubbles  are  shown  in  Fie.  h  for  >  -  2.3.  0.0.  and  0.42.  At 
v  =  0.42  the  opposite  sides  of  the  bubble  t ust  touched  one  another  on  the  v-axis  The 
bubble  might  then  split  into  two  bubbles  For  y  -  o  our  result  agreed  with  the  analytic 
solution  obtained  by  f  FI  Slel  eod  1 12,  in  this  special  ease. 

' '  'x 

v  -  U 


-0.2  il  I  n  ii  I  (1.2 


Fic ».  f>.  On,-  half  of  4  >„'<•,/  bubbU-  *->r  \  2  }.  \  -  u  otui  v  -  v..  --  -  n  4 2 

The  half-bubbles  shown  in  Fig.  b  also  represent  bubbles  attached  to  a  w  all  with 
contact  angle  0  =  rr' 2.  We  also  determined  the  shapes  of  bubbles  attached  to  a  wall  with 
other  values  of  0  by  a  suitable  modification  of  the  preceding  method  The  difference  in 
the  formulation  was  that  the  right  side  of  i5.4i  w  as  replaced  by  a:  tan  0.  Then  in  changing 
variables  we  set  -  1  -  n  ~  ‘ . 

For  each  value  of  0,  the  bubble  began  as  a  circular  are  of  radius  y  for  y  \  erv  lai  ge 
and  it  deformed  as  y  was  decreased.  Finally  a  value  of  y.  w  hich  ve  called  v.,igi  i,  w  as 
reached  at  which  opposite  sides  of  the  bubble  came  into  contact  w  ith  one  another.  For 
0  £  rr/2  the  contact  occurred  on  the  wall  as  in  Fig  6.  However  fo-  0  >  rr 1  2.  it  occurred 
off  the  wall  as  is  shown  in  Fig.  7  for  0  =  2  rr'3  =  120'.  Bubbles  are  presented  in  that 
figure  for  y  =  1.3,  y  =  -0.6  and  y  =  y,,1 2x-  '3 1  -  -  1 .6.  I  he  negative  values  of  >  indicate 
that  the  pressure  in  the  bubble  is  less  than  the  stagnation  pressure.  If  the  bubble  breaks 
into  two  at  the  point  of  contact,  part  of  the  bubble  will  be  off  the  wall  and  part  will 
remain  attached  to  the  wall,  and  our  calculation  determines  the  size  of  each. 
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Fig.  7.  Three  computed  profiles  for  a  bubble  on  o  wall  with  contact  angle  [I  =  1 20c  and  y  =  1 .3.  y  -  0.6 
and  y  -  y,,  =  - 1 .6.  At  y  =  y,>  =  -  l .6.  the  bubble  touches  itself  at  about  the  distance  v  =  0.35  from  the  wall. 


6.  Bubble  distortion  in  a  straining  flow.  In  the  mixing  of  two  fluids,  a  drop  of  one 
fluid  may  be  broken  up  into  smaller  drops  by  the  action  of  the  other  fluid  flowing  around 
it.  To  investigate  this  process  J.-M.  Vanden-Broeck  and  I  [  1 3]  considered  the  distortion 
of  a  drop  or  bubble  in  a  straining  flow,  a  problem  proposed  to  us  by  Professor  A. 
Acrivos.  As  in  the  preceding  section,  we  treated  the  two-dimensional  case  of  a  constant 
pressure  bubble  or  drop  in  an  in  viscid,  incompressible  fluid.  The  stream  function  of  the 
flow  was  assumed  to  be  ax y  far  from  the  bubble. 

The  formulation  of  this  problem  and  the  method  of  solution  were  very  similar  to 
those  of  the  preceding  section,  so  we  shall  just  indicate  the  differences  and  then  describe 
the  results.  The  length  scale  we  chose  was  (2tr/'pa‘)'  \  the  velocity  scale  was 
(2<Ta/p)'  \  and  then  the  sole  dimensionless  parameter  was  y  defined  by 

•  6-1)  y  =  2(  p,  -pb) Ultra  V  'p'  \ 

Because  of  the  symmetry  of  the  incident  flow,  it  sufficed  to  solve  the  problem  in  the 
angular  sector  0  S  <9  S  w/4  with  6  =  0  a  streamline  and  b  =  7t/4  a  potential  line.  We  then 
introduced  new  variables,  used  the  Hilbert  transform,  introduced  a  mesh  and  difference 
formulas,  etc.  Finally  we  solved  the  resulting  equations  by  Newton's  method. 

The  calculated  bubbles  for  three  values  of  y  are  shown  in  Fig.  8.  As  y-»oc,  the 
bubble  tends  to  a  circle  of  radius  y  and  the  bubble  shown  for  y  =  5  is  still  close  to  a 
circle.  However  when  y  has  decreased  to  zero,  the  bubble  is  nearly  a  square  with 
rounded  corners.  This  is  understandable  because  there  are  four  stagnation  points  where 
the  bubble  meets  the  x  and  v  axes.  At  them,  the  fluid  pressure  is  maximal  and  it  pushes 
the  bubble  surface  inward.  When  y  =  -1.24  these  four  points  have  moved  further  in 
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while  in  between  adjacent  pairs  of  them  four  horns  or  spikes  have  appeared.  We  have 
also  determined  the  ultimate  form  of  the  horns  as  y  is  further  decreased.  We  find  tha* 
near  each  tip  opposite  sides  come  into  contact  with  one  another  and  form  a  small  bubble 
there. 

7.  Damping  of  underwater  explosion  bubble  oscillations.  Now  we  shall  consider  a 
time  dependent  problem,  the  oscillation  of  the  gas  bubble  produced  by  the  detonation 
of  an  explosive  charge  under  water.  Such  a  bubble  is  initially  small,  spherical,  and  at 
very  high  pressure.  It  expands  rapidly,  remaining  spherical  and  pushing  the  water 
radially  outward.  W'hen  the  bubble  pressure  reaches  the  ambient  value,  the  water  is  still 
moving  outward,  causing  the  bubble  to  overexpand.  When  the  expansion  finally  stops, 
the  bubble  begins  to  contract  but  the  inertia  of  the  water  again  causes  it  to  overcontract 
and  the  process  repeats  itself.  In  this  way  the  radial  oscillations  of  underwater  explosion 
bubbles  are  produced.  The  dots  in  Fig.  9  are  experimental  values  of  the  radius  of  one 
such  bubble  at  different  times  after  detonation,  and  they  show  the  typical  oscillations. 

In  1916,  Rayleigh  formulated  the  theory  of  these  oscillations  on  the  assumption 
that  the  water  was  inviscid  and  incompressible.  His  theory  yielded  undamped  periodic 
oscillations  because  it  did  not  contain  any  mechanism  for  energy  loss.  It  was  later 
realized  that  acoustic  radiation  was  the  main  loss  mechanism.  This  was  verified  by 
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Fig.  9,  The  bubble  radius  as  a  function  of  time  fora  0.55  lb.  charge  of  tetryl  detonated  at  a  depth  of  300  ft. 
below  the  water  surface.  The  dots  are  experimental  points  and  the  solid  curve  is  calculated  from  (7.4).  The  unit 
of  length  is  S  =  6  inches  and  the  unit  of  time  is  dippo'Y'2  =  4.85  milliseconds. 
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estimating  the  acoustic  loss,  calculating  the  consequent  damping,  and  comparing  the 
result  with  observations. 

To  incorporate  the  loss  mechanism  into  the  theory,  I.  I.  Kolodner  and  1  [14],  in 
1 95 1 ,  developed  a  new  theory  of  the  oscillation  taking  account  of  the  compressibility  of 
the  water.  We  did  so  by  using  the  wave  equation  for  the  potential  function  of  the  water 
motion,  rather  than  Laplace's  equation  which  Rayleigh  had  used. 

The  use  of  the  wave  equation  to  describe  large  amplitude  motions  was  considered 
inappropriate  at  that  time,  and  it  still  is  so  considered  by  many  fluid  dynamicists.  On  the 
other  hand,  the  less  accurate  Laplace  equation  is  regarded  as  perfectly  suited  to 
describe  such  motions.  This  difference  in  attitude  is  probably  due  to  the  differences  in 
the  methods  of  derivation  of  the  two  equations.  In  the  derivation  of  the  wave  equation, 
the  smallness  of  the  velocity  is  emphasized  to  justify  linearization  of  the  compressibility 
terms.  However  in  the  derivation  of  Laplace  s  equation,  the  compressibility  terms  are 
just  ignored  from  the  beginning,  so  it  is  never  necessary  to  mention  that  the  velocity  is 
small. 

Our  theory  of  the  bubble  oscillation  is  based  on  the  assumption  that  the  pressure 
P(a )  in  the  bubble  is  a  function  of  the  bubble  radius  a  (r),  given  by  the  adiabatic  relation 

(7.1)  P(a)  =  /t[4fl’/3]  v. 

Here  k  and  y  are  given  constants.  The  potential  function  y>(r,  r)  of  the  water  motion  is 
assumed  to  satisfy  the  wave  equation  with  sound  speed  r. 

(7.21  A<p—  c  ~<f„  —  0. 

Then  the  pressure  p(r,  0  in  the  water  is  given  by  the  Bernoulli  equation 

(7.3)  p(r,  r)  =  p(,-p(<p,  +  yip,:). 

where  p0  is  the  initial  pressure  and  p  is  the  density  of  the  water.  The  formulation  is 
completed  by  requiring  the  pressure  p(a,  ft  in  the  fluid  to  equal  that  in  the  bubble,  P(a ), 
at  r  =  a,  and  requiring  the  fluid  velocity  <p,(a,  t)  to  equal  a,  at  r  =  a.  We  also  give  the 
initial  radius  and  velocity  of  the  bubble  and  assume  that  <p  =  <p,  =  0  for  r  >  a  at  i  =  0. 

From  this  formulation  we  derived  the  following  autonomous  nonlinear  second- 
order  ordinary  differential  equation  for  a(t): 

(7.4)  (a, -c)(aa„  +%a?  -A )-o?  +a  l(a2 A),  =  0. 

Here  A(a )  is  defined  by 

(7.5)  A(a)  =  p  l[f>(<i)-p„]. 

If  we  divide  (7.4)  by  c  and  let  c  become  infinite,  we  obtain  Rayleigh's  equation  for  a(/>. 
A  phase  plane  analysis  shows  that  the  equilibrium  point  of  ( 7.4 )  is  a  spiral  point, 
corresponding  to  damped  oscillations,  while  that  of  Rayleigh's  equation  is  a  center, 
representing  periodic  oscillations. 

In  Fig.  9  the  solid  curve  represents  the  solution  of  (7.4 1  obtained  numerically  for 
initial  conditions  corresponding  to  the  bubble  whose  observed  radius  is  given  by  the 
dots.  We  see  that  the  theoretical  curve  is  in  good  agreement  with  the  observed  data, 
and  it  correctly  predicts  the  damping  of  the  radial  motion  due  to  the  radiation  of 
acoustic  energy.  From  ait  I  the  pressure  pulses  radiated  by  the  bubble  can  be  found. 
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8.  Forced  oscillations  of  a  gas  bubble  in  a  sound  field.  To  describe  the  oscillation  of 
small  bubbles  in  a  fluid,  such  as  cav  itation  bubbles  or  air  bubbles,  it  is  necessary  to  take 
into  account  surface  tension  and  v  iscositv  as  well  as  inertia.  Therefore.  I’lesset  [15]  and 
others  modified  Rayleigh's  equation  to  include  these  erfect.,.  They  also  included  an 
external  time  dependent  pressure,  representing  an  incident  sound  field 

Lauterborn  [lb]  solved  the  modified  equation  numerically  for  periodic  incident 
sound  fields,  seeking  periodic  oscillations  of  the  bubble  He  obtained  them  for  small  and 
moderate  values  of  the  amplitude  of  the  applied  pressure.  From  his  results  he  plotted  a 
response  curve  for  each  applied  pressure,  showing  the  maximum  bubble  radius  of  the 
periodic  solution  as  a  function  of  the  applied  frequency.  However,  for  very  large 
amplitudes  of  the  applied  pressure,  the  solution  did  not  become  periodic;  or  if  it  did,  the 
amplitudes  did  not  lie  on  a  smooth  curve. 

It  seemed  to  me  that  this  difficulty  with  the  large  amplitude  oscillations  could  be 
avoided  by  taking  into  account  the  acoustic  radiation  from  the  bubble,  just  as  was  done 
in  the  preceding  section  to  describe  the  damping  of  underwater  explosion  bubble 
oscillations.  Therefore  my  student  M.  Miksis  and  I  [17]  modified  1 7.4 1,  the  equation  of 
Keller  and  Kolodner,  to  include  surface  tension,  viscosity  and  an  incident  sound  field,  in 
the  same  way  as  Plesset  and  others  had  modified  Rayleigh's  equation. 

The  equation  we  derived  in  this  way  is 


(8.1i 
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Here  ait)  is  the  bubble  radius  at  time  t ,  c  is  the  sound  speed  in  the  fluid,  p  is  the  fluid 
density,  p  is  the  coefficient  of  v  iscosity  of  the  fluid,  <r  is  the  surface  tension,  2c  g '  is  the 
time  derivative  of  the  potential  of  the  incident  sound  field  evaluated  at  the  center  of  the 
bubble  and  pMa)  is  the  difference  between  the  bubble  pressure  and  the  pressure  at 
infinity.  The  pressure  in  the  bubble  was  taken  to  be  ka  ,v  *  pt,  where  p„  represents  a 
constant  vapor  pressure  and  k  and  y  are  constants.  When  c  =  x,  (8.1 )  reduces  to  the 
equation  solved  by  Lauterborn;  while  when  p  =  v  =  g  -  p,  =0,  it  reduces  to  the 
equation  (7.4)  of  Keller  and  Kolodner. 

To  treat  this  equation  we  first  studied  the  free  oscillations  corresponding  to  g  -  o 
by  a  phase  plane  analysis,  and  obtained  results  like  those  of  Keller  and  Kolodner  [14], 
We  also  solved  (8.11  numerically  for  the  same  case  shown  in  Fig.  9  and  obtained  the 
same  result  because  surface  tension  and  viscosity  are  negligible  in  that  case. 

Next  we  chose  2c  V'(rt  =  p  'f*  sin  ml,  which  is  the  same  forcing  function  as  was 
used  by  Lauterborn  [16],  Then  we  solved  the  equation  analytically  by  the  method  of 
averaging,  in  the  same  way  as  Prosperetti  [18]  had  solved  the  equation  used  bv 
Lauterborn.  The  results  were  similar  to  his,  showing  harmonic,  subharmonic  and 
ultraharmonic  resonances. 

Finally  we  solved  the  equation  numerically  as  an  initial  value  problem,  as  Lauter¬ 
born  did.  Some  of  our  results  are  shown  in  Fig.  10.  They  agree  with  his  except  at  the 
largest  forcing  amplitudes.  Then  our  equation  yields  periodic  solutions  and  a  response 
curve  similar  to  those  which  we  both  found  for  smaller  amplitudes,  whereas  his  did  not. 
Thus  the  inclusion  of  the  effects  of  acoustic  radiation  in  (8. 1  >  has  had  the  desired  effect 
of  yielding  periodic  oscillations  and  a  response  curve  for  even  very  large  amplitudes  of 
the  incident  pressure  field. 
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